% I0是输入图像，I0_name是输出图像名，rgbia(1为r,2为g，3为b),spd是物像素，spdh是全息图像素，z是衍射距离
function dot2plane(I0,I0_name,rgbia,rgbim,rgbib,spd,spdh,z,rh,ch)

   for d2pi = rgbia:rgbim:rgbib
    
%======================================================================================
%数据处理

I1 = im2double(I0(:,:,d2pi)); 
% I1=im2double(I0(1:4,1:4,1));        %取一层的值，如r、g、b之一  
lamb_3 = [0.000650 0.000520 0.000450];
lamb_I = lamb_3(d2pi);


[r,c] = size(I1); 
Amp = zeros(r*c,1);
xyDis = zeros(r*c,2);
I1_rand = I1.*exp(1j*rand(r,c)*2*pi);

 I2 = zeros(r,c);

% spd = 0.1; %物面抽样距离,
spd = spd;

ymo = mod(r,2);
xmo = mod(c,2);
num = 1;
for i=1:r
    yod = -(r-ymo)/2*spd+i*spd;
    for n=1:c
        if I1(i,n) > 0
            I2(i,n) = I1(i,n);
            xod = -(c-xmo)/2*spd+n*spd;            
            Amp(num) =  I1_rand(i,n);%光强
            xyDis(num,1) = xod;%x位置
            xyDis(num,2) = yod;%y位置
            num = num+1;
        end
    end
end

o_amp = Amp(1:num-1,:);
o_xydis = xyDis(1:num-1,:);


%===============================================================================
%全息图计算

[o_num,amp_num] = size(o_amp);

lamb = lamb_I;        %激光波长，单位mm
k = 2 * pi / lamb;
% z = 300;
z = z;
% z2_1 = 1/2/z;

%====全息图参数=============================================

%全息图点距
% spdh = 3.18e-4;
spdh = spdh;

%全息图点数
% rh = 10000;
% ch = 10000;
rh = rh;
ch = ch;

 %全息图大小
Lx=spdh*ch;  
Ly=spdh*rh;

%全息图偏置量
yh_shift3 = [-58 -47 -35]; 
xh_shift = 0; 
yh_shift = yh_shift3(d2pi); 

xh = linspace(-Lx/2,Lx/2,ch);
yh = linspace(-Ly/2 + yh_shift,Ly/2 + yh_shift,rh);

[xh,yh] = meshgrid(xh,yh);


jkz = 1j*k*z;
jlz = 1j*lamb*z;
expjkz_jlz = exp(jkz)/jlz;

jk_2z = 1j*k/2/z;

UF = zeros(rh,ch);
UF1 = zeros(rh,ch);

   for onum = 1:o_num
      tic
    
      xo = o_xydis(onum,1);
      yo = o_xydis(onum,2);

%发散球面光波
% UF1 = exp(1j*k*z).*exp(1j*k.*((xh-xo).^2+(yh-yo).^2)/2/z);  %发散球面光波
% UF1 = exp(jk.*xyz)./xyz;
    UF1 = o_amp(onum,1)*expjkz_jlz.*exp(jk_2z.*((xh-xo).^2+(yh-yo).^2));  %发散球面光波

%会聚球面光波
% UF1 = exp(-1j*k*z).*exp(-1j*k.*((xh-xo).^2+(yh-yo).^2)/2/z);  %会聚球面光波

    UF = UF + UF1;

    toc
    disp(onum/o_num);
   end
   
save("UF"+I0_name+d2pi+".mat","UF");

%参考光
xR = linspace(-Lx/2,Lx/2,ch);
yR = linspace(-Ly/2,Ly/2,rh);
[xR,yR] = meshgrid(xR,yR);

R=k*yR*sind(0)+k*xR*sind(0); 
II=abs(UF).*cos(angle(UF)-R);
II(II>=0)=255;
II(II<0)=0;

II=(II-min(min(II)))./(max(max(II))-min(min(II)))*255;
[AAA,map]=imread("tg1-gray.bmp");
imwrite(II,map,I0_name+d2pi+".bmp");

   end
end


